function [q_T, q_T_F, E_tax_F, E_tax_T_F, E_MC_tax_F, E_MC_tax_T_F, max_MC_tax] = ...
    fn_E_tax(R1,delta,D0,chi,phi,s_min,~,dist_s,dist_d,d_min,d_max,prob,s_hat,s_star,kappa)

%% Probability of failure

q_F   = dist_s.cdf(s_hat) + prob*(dist_s.cdf(s_star) - dist_s.cdf(s_hat));

%% Probability of paying tax (conditional on failure)

f_q_T_F = @(s) (fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max)>0).*dist_s.pdf(s);

T1 = integral(f_q_T_F,s_min,s_hat);
T2 = integral(f_q_T_F,s_hat,s_star);

q_T   = T1 + prob*T2;
q_T_F = q_T/q_F;

%% Expected tax, conditional on failure, and conditional on failure and T>0

f_T_F = @(s) fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max).*dist_s.pdf(s);

T1 = integral(f_T_F,s_min,s_hat);
T2 = integral(f_T_F,s_hat,s_star);

E_tax_F   = (T1 + prob*T2)/q_F;
E_tax_T_F = (T1 + prob*T2)/q_T;

%% Expected marginal cost of tax, conditional on failure, and conditional on failure and T>0

f_MC_F = @(s) kappa(1)*exp(kappa(2)*fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max)).*(fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max)>0).*dist_s.pdf(s);

T1 = integral(f_MC_F,s_min,s_hat);
T2 = integral(f_MC_F,s_hat,s_star);

E_MC_tax_F   = (T1 + prob*T2)/q_F;
E_MC_tax_T_F = (T1 + prob*T2)/q_T;

max_MC_tax = kappa(1)*exp(kappa(2)*fn_tax(delta,R1,D0,chi,s_min,phi,dist_d,d_min,d_max));

% Note: here we compute the expectation of the marginal cost, conditional
% on failure. We could compute instead the marginal cost, conditional on
% failure and T>0. This  would imply that the mg. cost in calibration_norm
% becomes: q_F_1*q_F_T_1*E_MC_tax_pos_1*unins_acc_1

end